home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / int_tabulated.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  172 lines
  1. ; $Id: int_tabulated.pro,v 1.13 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1995-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       INT_TABULATED
  8. ;
  9. ; PURPOSE:
  10. ;       This function integrates a tabulated set of data { x(i) , f(i) },
  11. ;       on the closed interval [min(X) , max(X)].
  12. ;
  13. ; CATEGORY:
  14. ;       Numerical Analysis.
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = INT_TABULATED(X, F)
  18. ;
  19. ; INPUTS:
  20. ;       X:  The tabulated X-value data. This data may be irregularly
  21. ;           gridded and in random order. If the data is randomly ordered
  22. ;        you must set the SORT keyword to a nonzero value.
  23. ;           Duplicate x values will result in a warning message.
  24. ;       F:  The tabulated F-value data. Upon input to the function
  25. ;           X(i) and F(i) must have corresponding indices for all
  26. ;        values of i. If X is reordered, F is also reordered.
  27. ;
  28. ;       X and F must be of floating point or double precision type.
  29. ;
  30. ; KEYWORD PARAMETERS:
  31. ;       SORT:   A zero or non-zero scalar value.
  32. ;               SORT = 0 (the default) The tabulated x-value data is
  33. ;                        already in ascending order.
  34. ;               SORT = 1 The tabulated x-value data is in random order
  35. ;                        and requires sorting into ascending order. Both
  36. ;             input parameters X and F are returned sorted.
  37. ;       DOUBLE: If set to a non-zero value, computations are done in
  38. ;               double precision arithmetic.
  39. ;
  40. ; OUTPUTS:
  41. ;       This fuction returns the integral of F computed from the tabulated
  42. ;    data in the closed interval [min(X) , max(X)].
  43. ;
  44. ; RESTRICTIONS:
  45. ;       Data that is highly oscillatory requires a sufficient number
  46. ;       of samples for an accurate integral approximation.
  47. ;
  48. ; PROCEDURE:
  49. ;       INT_TABULATED.PRO constructs a regularly gridded x-axis with a
  50. ;    number of segments as an integer multiple of four. Segments
  51. ;    are processed in groups of four using a 5-point Newton-Cotes
  52. ;    integration formula.
  53. ;       For 'sufficiently sampled' data, this algorithm is highly accurate.
  54. ;
  55. ; EXAMPLES:
  56. ;       Example 1:
  57. ;       Define 11 x-values on the closed interval [0.0 , 0.8].
  58. ;         x = [0.0, .12, .22, .32, .36, .40, .44, .54, .64, .70, .80]
  59. ;
  60. ;       Define 11 f-values corresponding to x(i).
  61. ;         f = [0.200000, 1.30973, 1.30524, 1.74339, 2.07490, 2.45600, $
  62. ;              2.84299,  3.50730, 3.18194, 2.36302, 0.231964]
  63. ;
  64. ;       Compute the integral.
  65. ;         result = INT_TABULATED(x, f)
  66. ;
  67. ;       In this example, the f-values are generated from a known function,
  68. ;       (f = .2 + 25*x - 200*x^2 + 675*x^3 - 900*x^4 + 400*x^5)
  69. ;
  70. ;       The Multiple Application Trapazoid Method yields;  result = 1.5648
  71. ;       The Multiple Application Simpson's Method yields;  result = 1.6036
  72. ;                INT_TABULATED.PRO yields;  result = 1.6232
  73. ;         The Exact Solution (4 decimal accuracy) yields;  result = 1.6405
  74. ;
  75. ;    Example 2: 
  76. ;       Create 30 random points in the closed interval [-2 , 1].
  77. ;         x = randomu(seed, 30) * 3.0 - 2.0
  78. ;
  79. ;       Explicitly define the interval's endpoints.
  80. ;         x(0) = -2.0  &  x(29) = 1.0
  81. ;
  82. ;       Generate f(i) corresponding to x(i) from a given function.
  83. ;         f = sin(2*x) * exp(cos(2*x))
  84. ;
  85. ;       Call INT_TABULATED with the SORT keyword.
  86. ;         result = INT_TABULATED(x, f, /sort)
  87. ;
  88. ;       In this example, the f-values are generated from the function,
  89. ;       f = sin(2*x) * exp(cos(2*x))
  90. ;
  91. ;       The result of this example will vary because the x(i) are random.
  92. ;       Executing this example three times gave the following results:
  93. ;                       INT_TABULATED.PRO yields;  result = -0.0702
  94. ;                       INT_TABULATED.PRO yields;  result = -0.0731
  95. ;                       INT_TABULATED.PRO yields;  result = -0.0698
  96. ;        The Exact Solution (4 decimal accuracy) yields;  result = -0.0697
  97. ;
  98. ; MODIFICATION HISTORY:
  99. ;           Written by:  GGS, RSI, September 1993
  100. ;           Modified:    GGS, RSI, November  1993
  101. ;                        Use Numerical Recipes cubic spline interpolation 
  102. ;                        function NR_SPLINE/NR_SPLINT. Execution time is 
  103. ;                        greatly reduced. Added DOUBLE keyword. The 'sigma' 
  104. ;                        keyword is no longer supported.
  105. ;           Modified:    GGS, RSI, April  1995
  106. ;                        Changed cubic spline calls from NR_SPLINE/NR_SPLINT
  107. ;                        to SPL_INIT/SPL_INTERP. Improved double-precision
  108. ;                        accuracy.
  109. ;           Modified:    GGS, RSI, April 1996
  110. ;                        Replaced WHILE loop with vector operations. 
  111. ;                        Check for duplicate points in x vector.  
  112. ;                        Modified keyword checking and use of double precision.
  113. ;-
  114.  
  115. FUNCTION Int_Tabulated, X, F, Double = Double, Sort = Sort
  116.  
  117.   ;Return to caller if an error occurs.
  118.   ON_ERROR, 2 
  119.  
  120.   TypeX = SIZE(X)
  121.   TypeF = SIZE(F)
  122.  
  123.   ;Check F data type.
  124.   if TypeF[TypeF[0]+1] ne 4 and TypeF[TypeF[0]+1] ne 5 then $
  125.     MESSAGE, "F values must be float or double."
  126.  
  127.   ;Check length.
  128.   if TypeX[TypeX[0]+2] ne TypeF[TypeF[0]+2] then $
  129.     MESSAGE, "X and F arrays must have the same number of elements."
  130.  
  131.   ;Check duplicate values.
  132.   if TypeX[TypeX[0]+2] ne N_ELEMENTS(UNIQ(X[SORT(X)])) then $
  133.     MESSAGE, "X array contains duplicate points."
  134.  
  135.   ;If the DOUBLE keyword is not set then the internal precision and
  136.   ;result are identical to the type of input.
  137.   if N_ELEMENTS(Double) eq 0 then $
  138.     Double = (TypeX[TypeX[0]+1] eq 5 or TypeF[TypeF[0]+1] eq 5) 
  139.  
  140.   Xsegments = TypeX[TypeX[0]+2] - 1L
  141.  
  142.   ;Sort vectors into ascending order.
  143.   if KEYWORD_SET(Sort) ne 0 then begin
  144.     ii = SORT(x)
  145.     X = X[ii]
  146.     F = F[ii]
  147.   endif
  148.  
  149.   while (Xsegments MOD 4L) ne 0L do $
  150.     Xsegments = Xsegments + 1L
  151.  
  152.   Xmin = MIN(X)
  153.   Xmax = MAX(X)
  154.  
  155.   ;Uniform step size.
  156.     h = (Xmax+0.0 - Xmin) / Xsegments
  157.   ;Compute the interpolates at Xgrid.
  158.     ;x values of interpolates >> Xgrid = h * FINDGEN(Xsegments + 1L) + Xmin
  159.     z = SPL_INTERP(X, F, SPL_INIT(X, F, Double = Double), $
  160.                    h * FINDGEN(Xsegments + 1L) + Xmin, Double = Double)
  161.   ;Compute the integral using the 5-point Newton-Cotes formula.
  162.     ii = (LINDGEN((N_ELEMENTS(z) - 1L)/4L)+1) * 4
  163.     if Double eq 0 then $
  164.       RETURN, FLOAT(TOTAL(2.0 * h * (7.0 * (z[ii-4] + z[ii]) + $
  165.                     32.0 * (z[ii-3] + z[ii-1]) + 12.0 * z[ii-2]) / 45.0)) $
  166.     else $
  167.       RETURN, TOTAL(2D * h * (7D * (z[ii-4] + z[ii]) + $
  168.                     32D * (z[ii-3] + z[ii-1]) + 12D * z[ii-2]) / 45D, /DOUBLE)
  169.  
  170. END
  171.  
  172.